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To develop a full understanding of interactions in nanomagnet arrays is a persistent challenge, critically 
impacting their technological acceptance. This paper reports the experimental, numerical and analytical 
investigation of interactions in arrays of Co nanoellipses using the first-order reversal curve (FORC) 
technique. A mean-field analysis has revealed the physical mechanisms giving rise to all of the observed 
features: a shift of the non-interacting FORC-ridge at the low-H c end off the local coercivity H c axis; a 
stretch of the FORC-ridge at the high-if c end without shifting it off the H c axis; and a formation of a 
tilted edge connected to the ridge at the low-H c end. Changing from fiat to Gaussian coercivity 
distribution produces a negative feature, bends the ridge, and broadens the edge. Finally, nearest 
neighbor interactions segment the FORC-ridge. These results demonstrate that the FORC approach 
provides a comprehensive framework to qualitatively and quantitatively decode interactions in 
nanomagnet arrays. 

Nanomagnet arrays are basic building blocks 1 for key technologies such as ultrahigh density magnetic 
recording media 2 " 4 , magnetic random access memory (MRAM) 5,6 , and logic devices 7 " 9 . Interactions 
within the arrays critically affect the functionalities of the nanomagnets as well as enable new device 
concepts. For instance, dipolar interactions may trigger an analog memory effect in nanowire arrays 10 , enable 
digital computation in magnetic quantum-dot cellular automata systems 7,811 , lead to frustrations in artificial "spin 
ice" 12 " 14 , or adversely affect thermal stability and switching field distribution in magnetic recording media or 
MRAM elements 5,15 " 17 . Probing and managing these interactions is often difficult because they are long-ranged, 
anisotropic, and configuration-dependent 17 . 

The first-order reversal curve (FORC) method 18,19 has provided detailed characterization for a variety of 
magnetic 20 " 27 and other hysteretic systems 28,29 . However, a coherent framework to interpret the features of the 
FORC diagrams and extract quantitative information is still lacking despite decades of effort by numerous groups. 
In this work, using the FORC method we have quantitatively investigated tunable interactions in model systems 
of single domain nanomagnet arrays. With mean-field level simulations, supplemented with a cluster extension, 
we have reproduced all features and trends of the experimental FORC diagrams quantitatively and identified their 
physical origins. Our approach decodes interactions in nanomagnet arrays, even disordered arrays, and also 
presents a pathway to evaluate (de) stabilizing interactions in other hysteretic systems. 

Results 

Rectangular arrays of polycrystalline Co ellipses were fabricated with varying center-to-center spacing by mag- 
netron sputtering, in conjunction with e-beam lithography and lift-off. Details are presented in Methods. The 
ellipses have major/minor axes of 220/1 10 nm, with a structure of Ta(l nm)/Co(9 nm)/Ta(l nm), forming 50 X 
50 um 2 arrays. In arrays Al/2/3 the minor-axis spacings of 150/200/250 nm are less than the major-axis spacing 
of 500 nm. Thus the mean dipolar interactions are demagnetizing, favoring anti-parallel alignment. In arrays Bl/ 
2/3 the minor-axis spacing of 500 nm exceeds the major-axis spacings of 250/300/350 nm. Therefore, the mean 
dipolar interactions are magnetizing, favoring parallel alignment. Scanning electron microscopy (SEM) and 
magnetic force microscopy (MFM) images, at remanence after DC demagnetization, of arrays Al and Bl are 
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Figure 1 | Demagnetizing Arrays: (a) SEM and (b) MFM image of the DC- 
demagnetized Al array. Dashed ovals outline single ellipses, while the 
dashed box outlines an example of the checkerboard pattern, (c) 
Experimental and (d) simulated FORC distributions for the A2 array. 
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Figure 2 | Magnetizing Arrays: (a) SEM and (b) MFM image of the 

DC-demagnetized Bl array. Dashed ovals outline single ellipses. 

(c) Experimental and (d) simulated FORC distributions for the B2 array. 



shown in Figs. 1 and 2, respectively. MFM image contrast indicates 
the out-of-plane stray fields and confirms the ellipses' single domain 
state. 

FORC measurements were performed to obtain magnetization 

M(H, H R ) under different reversal field H R and applied field 

H. The FORC-distribution is then extracted 18 , p(H,H R ) = 

1 d 2 M(H,H R ) w . t 
— , where M s is the saturation magnetization. 

2M S dHdH R 5 

In loose analogy to the Preisach model, the FORC-distribution in 
certain simple cases can be interpreted as the 2-dimensional dis- 
tribution of elemental hysteresis loops with unit magnetization 
called hysterons on the (H, H R ) plane, or on the corresponding 
(H c , H B ) plane, defined by local coercivity H c = (H-H R )/2 and 
bias/interaction field H B = (H+H R )/2. 

Good agreement between measured and simulated FORC- distri- 
butions was obtained for all studied arrays [e.g., Figs. 1(c) and (d) for 
array A2]. All demagnetizing FORCs A 1/2/3 exhibited a ridge with 
the high-H c end on the H c axis and the low-H c end shifted in the 
+H B direction 30 . Increasing interactions in A3^A2^>A1 increased 
the low-H c shift and the length of the ridge. In addition, an edge 
emerged from the low-H c end towards negative H B , highlighted by 
the arrows, forming a "wishbone" or boomerang structure 24 . A nega- 
tive feature is observed at negative H B values near the high-H c end. 
Similarly, FORC distributions for the magnetizing arrays also exhibit 
a ridge with the high-H c end on the H c axis, but with the low-H c end 
shifted towards — H B [e.g., Figs. 2(c) and (d) for array B2]. Increasing 
interactions again increase the low-H c shift, but reduce the length of 
the ridge 30 . A negative feature below the ridge is more prominent 
than that in the demagnetizing case. 

Non-interacting case. A FORC gives a non-zero contribution to 
p(H, H R ) only if dM/dH along the FORC depends on H R . We first 
show that in the non-interacting case p coincides with the coercivity 
distribution D(H K ), spread along the H c axis. Indeed, particle P z with 
coercivity H K * down-flips at H dn { = —H K ' and up-flips at H u p = H K \ 
Therefore, on a FORC starting at H R > —H K \ P { starts out up-flipped 
and remains so, not contributing to dM/dH nor p. 

In contrast, on FORC(H, H R = —Hk), Pi is the last to down-flip 
along the major loop, and has the highest coercivity among the 
down-flipped particles. Therefore, P { is the last to up-flip as H 
increases past H K { on the same FORC, causing a dM/dH > 0 jump. 
Since this dM/dH jump is unmatched by the neighboring FORC(H, 
H R > —H^), dM/dH exhibits a dependence on H R , making d(dM/ 



dH)/dH R non-zero. dM/dH increases as H R decreases, making a pos- 
itive contribution to p at (H = H K \ H R = —Hk). 

For all subsequent FORCs at H R < —H K \ P { starts out down- 
flipped but still up-flips at H = H K \ The dM/dH jumps on these 
FORCs are matched since they occur at the same field on each 
FORC(H, H R < -H^). Thus dM/dH is independent of H R , and 
doesn't contribute to p. This reasoning highlights that only dM/dH 
jumps on individual FORCs that are unmatched by neighboring 
FORCs contribute to p. Each particle P { contributes to p only once, 
at (H = H K \ H R = -Hk 1 ) or equivalently at (H c = H K \ H B = 0). The 
contributions of all particles gives rise to a ridge along the H c (H R = 
—H) axis, which reflects D(H K ). 

Interacting case. Next, we introduce interactions between nano- 
magnets on the mean-field level by including an interaction field 
H int = ocM(H), where a<0 for demagnetizing systems and a>0 for 
magnetizing ones 21 . Fig. 3(a) shows a sequence of FORCs for a 
demagnetizing system with a rectangular D(H K ), and a zoom-in 
view of the boxed region (right panel). The three FORC segment- 
pairs (l)/(2)/(3) show that the last dM/dH jump on each FORC(H, 
H R = —HK-aMiHp)) - caused by the last up-flipping particle P(Hk) 
- is unmatched by the neighboring FORC(H, H R > — Hk-olM(Hr)). 

Fig. 3(b) shows that with interactions the unmatched dM/dH 
jumps still generate the ridge, but at shifted H B values. 
Importantly, on the mean-field level all particles experience the same 
interaction field and thus the order of flips continues to be governed 
by the order of the coercivities: along the major loop the particles 
down-flip in the order of their coercivities, lowest (highest) coercivity 
particle P(H K min ) first [P(H K max ) last] . Starting at the low-H c end, the 
lowest coercivity particle P(H K min ) down-flips at H dn min = —H K min - 
aM s , as no other particles have flipped yet: M(H dn min ) = M s . 
Increasing H along FORC(H, H R = H dn min \ P(H K min ) up-flips at 
H up min = H K min -otM s , causing a positive jump dM/dH > 0 as shown 
by the lower FORC(H, H R = H dn min ) -segment of pair (1). This jump, 
caused by P{H K mi % is absent on the upper FORC(H, H R > H dn min )- 
segment and is thus unmatched, contributing to p at (H = H K min - 
aM s , H R = -H K min -OLM s ), or similarly at (H c = H K min , H B = -oeM s ), 
defining the low-H c end. These flipping fields are shifted from their 
non-interacting values, as shown by the arrow set (1) in Fig. 3(b). 
Since P(H K min ) defines the low-H c end of the FORC-ridge, one con- 
cludes that interactions shift the low-H c end of the FORC-ridge to 
the +H B direction by — otM s (recall a < 0), but leave its H c coord- 
inate un-shifted at H c = H K min . 
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Figure 3 | (a) Schematic illustration of family of FORCs for arrays with a flat coercivity distribution and mean-field demagnetizing interactions, 
with bold lines and numbers indicating unmatched dM/dH jumps. Calculated FORC distributions are shown in (b) illustrating the construction of the 
ridge and (c) the edge, (d) FORC distribution with the same interactions, but a Gaussian coercivity distribution; emergent negative feature is indicated by 
the dashed boundary. Similar panels are shown in (e-h) for the magnetizing case. 



FORC-segment pairs (2)/(3) in Fig. 3(a) illustrate the up-flip of 
higher coercivity particles P{H^). P{H^)s create unmatched dM/dH 
jumps on FORCs where they were the last to down-flip at H dn = 
—HK-ocMfHn) [vertical arrows in Fig. 3(b)] and also are the last to 
up-flip at H up = Hk -otM s [horizontal arrows in Fig. 3(b)], since 
M(H up ') = M s when P(h/) up-flips. 

The high-H c end of the FORC-ridge is defined by P(H K max ) that is 
the last to down-flip when the rest of the system is already negatively 
saturated (M = -M s ), thus H dn max = -H K max -(-otM s ). P(H K max ) up- 
flips along the FORC(H, H R = H dn max ) only after the rest of the 
system is positively saturated: H up max = H K max -ocM s . Accordingly, 
the unmatched dM/dH jumps caused by P(H K max ) contributes to p 
only at (H = H K max -otM s , H R = -H K max + oM s ), or similarly at (H c = 
H K max -aM s , H B = 0), defining the high-H c end. As observed before, 
the high-H c end of the FORC-ridge remains on the H c axis 24 , but 
stretched along the H c axis by — ocM s (a<0). 

Note that interactions shift the FORC ridge feature in H uniformly 
[Fig. 3(b)], i.e., the resultant projection of FORC distribution onto 
the H-axis is only displaced from its intrinsic values, but not dis- 
torted. This H-projection therefore mirrors the non-interacting case, 
where H up = H c = H, reflecting the intrinsic coercivity distribution, 
simply displaced by — otM s . Thus the intrinsic coercivity distribution 
can be - without distortion from interactions - directly identified 
from the FORC distribution. 

Figs. l(c,d) and 3(b,c) show that besides the ridge, p exhibit an 
edge as well with interactions 27 . As discussed earlier, in the absence of 
interactions, the dM/dH jumps along a FORC(H, H R = —H^) are 
matched by the jumps on the subsequent FORC(H, H R < —Hk)s, 
not contributing to p. The arrows of Fig. 3(a) show that the interac- 
tions destroy this matching specifically at the low-H c end by shifting 
the first up-flip field H up min of each FORC, caused by P{H K min ) y by 
— otMfHn). These shifts make the dM/dH jumps misaligned, see 
FORC-segment-pairs (4) and (5) (above and below), thus contrib- 
uting to p at (H = H^-aMiHp), H R ). These unmatched jumps give 
rise to the edge in Fig. 3(c). The end-points of the edge are (H = 
H K min -otM s , H R = —H K min -otM s ) and (H = H K min +aM s , H R = 
-H K max +otM s ), or alternatively (H c = H K min , H B = -oM s ) 
and (H c = (H K min +H K max )/2, H B = ocM s + (H K min -H K max )/ 2) . 
Accordingly, the tilt and asymmetry of the edge provide a direct 
measure of the width of coercivity distribution D(H K ) = H K max - 



H K min . In the extreme case of nearly identical nanomagnets with 
D(H K ) = H K max -H K min ^0, the edge is vertical in the H c -H B plane, 
at H c = H K and H B within ± aM s , as observed experimentally in Ni 
nanowire arrays 10 . 

In short, on the mean-field level the FORC-distribution of a sys- 
tem with demagnetizing interactions exhibits an edge and a ridge, 
shifted by the unmatched first and last dM/dH jumps along each 
FORC. The FORC-distribution vanishes between them for the con- 
sidered flat coercivity distribution, because jumps between the first 
and the last jumps along each FORC are matched by jumps on the 
neighboring FORCs. Here, the matched jumps are not caused by the 
same particles, as in the non-interacting system, but rather by dif- 
ferent particles whose up -flipping fields were shifted into alignment 
by the interactions. Still, the jumps are matched because the values of 
the aligned jumps are the same on neighboring FORCs for a flat 
distribution D(H K ). Visibly, the flat coercivity distribution on the 
mean-field already reproduces most features of the measured 
FORC-distribution. 

To improve our model we introduce a more realistic Gaussian 
D(H K ) to elucidate the origin of the negative features, which repres- 
ent a clear distinction between FORC and a literal Preisach inter- 
pretation. A Gaussian breaks the matching of jumps as now shown 
by examining the set of particles {P(H K Cent )} around the center of the 
coercivity distribution. A P(H K Cent ) particle is the last to down-flip at 
H R = -H K Cen \ where MfHj = 0, and the last to up-flip on the 
FORC(H, H R = -H K Cent ), contributing to p at (H = H K Cent -otM s , 
H R = -H K Cent ). On ¥ORC(H, H R < -H K Cent )s P{H K Cent ) is no longer 
the last to up-flip. On subsequent FORCs, up-flip jumps from dif- 
ferent particles get shifted into alignment with this jump. For the flat 
distribution, the number of particles shifted into alignment is steady, 
making the match complete and thus zero contribution to p. 
However, for the Gaussian distribution, the particles shifted into 
alignment come from the decreasing slope of the Gaussian, leading 
to dM/dH jumps with a decreasing magnitude, providing only a 
partial-match and generating a negative contribution to p. Fig. 3(d) 
shows the FORC-distribution for a Gaussian model that indeed 
develops a negative region specifically tracking the decreasing slope 
of the Gaussian, highlighted by the dashed line. Analogous argu- 
ments show that a Gaussian D(H K ) also bends the ridge and broadens 
the edge. 
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Figure 4 | Calculated (open symbols) and experimentally determined (solid symbols) interaction field for (a) demagnetizing arrays A 1/2/3 and (b) 
magnetizing arrays Bl/2/3. Averaged FORC- distribution utilizing the H c (threshold) are shown in insets for the (a) A2 and (b) B2 array, where the linear 
extrapolation is illustrated by the dashed line and open circle. 



Nearest-neighbor correlations. The last unaccounted feature of the 
measured p is the segmenting of the FORC-ridge into separate low- 
He and high-He ends, with different amplitudes. To explain this we 
include the nearest-neighbor interaction fields H nn l (conf) as the first 
terms of a systematic cluster-expansion. 

Decreasing the field from positive saturation, the weakest coercivity 
particles down-flip first. For demagnetizing interactions, H nr l(conf) of 
these down-flipped particles stabilize their nearest-neighbors in their 
up state. Therefore, for a sufficiently narrow D(Hk), the magnetiza- 
tion decreases towards zero by developing a checkerboard pattern 
[Fig. 1(b)]. The checkerboard naturally forms defects where the 
sequence of increasing coercivities selects third- nearest neighbor part- 
icles to down-flip. Still, the dominant reversal mechanism for nearly 
half of the particles is the checkerboard formation: down-flipping with 
all-up neighbors. Accordingly, the (nearly) half of the FORC ridge 
with Hk < H K Cent gets shifted along the +H B axis by H nr l(up), where 
Hnn(up) is the interaction field for the all-neighbors-up configuration. 
Once the checkerboard pattern is formed, the rest of the particles flip 
with neighbors in various intermediate configurations. Therefore, the 
Hk > H K Cent half of the ridge is broken into several pieces, shifted by 
varying H nr i(conj) fields. Consequently, the nearest neighbor interac- 
tions manifest as segmenting of the ridge 30 . 

To reiterate, the demagnetizing interactions (a < 0) (1) shift the 
non-interacting FORC-ridge at the low-H c end to the +H B direction 
by — otM s ; (2) stretch the non-interacting FORC-ridge at the high-H c 
end along H c by — ocM s without shifting it off the H c axis; and (3) 
form a tilted edge connected to the ridge at the low-H c end. 
Changing from flat to Gaussian D(Hk) distribution (4) produces a 
negative feature, bends the ridge, and broadens the edge. Finally, (5) 
nearest neighbor interactions segment the FORC-ridge. 

Magnetizing interactions. Adapting the above arguments for 
magnetizing interactions (a > 0): (1) the low-H c end is shifted in 
the — H B direction, and (2) the high-H c end is compressed without 
shifting it off the H c axis [Fig. 3(f)]. (3) Regarding the edge, the first 
up-flip along each branch is shifted by interactions in the opposite 
direction as the demagnetizing case [Fig. 3(e) right panel] . Therefore, 
the first dM/dH jumps are unmatched, decreasing in magnitude with 
more negative H R , thus negatively incrementing the FORC, forming a 
negative edge [Fig. 3(g)]. Changing from flat to Gaussian D(Hk) 
distribution (4) the negative edge gets pressed towards the positive 
ridge, and the FORC-ridge becomes curved [Fig. 3(h)]. The inclusion 
of nearest neighbor terms leads to (5) an avalanche reversal, 
collapsing the FORC-ridge to a single-value 31 . 

Quantifying Interaction Fields. Finally, we demonstrate the quan- 
titative predicting power of the above considerations. The lowest H K \ 
which is shifted in H B by aM s , is extracted from the FORC ridge by 
selecting an H c (threshold) such that 10% of the particles have H K l < 
H c (threshold), and averaging p over the H c = 0^H c (threshold) 
range. The averaged (dM/dH B ) y are shown in insets of Figs. 4(a) 
and 4(b). The H B shift is determined by linearly extrapolating 



(dM/dH B )' at the high \H B \ end to zero. The interaction field is 
calculated by a finite element method (using the NIST OOMMF 
code) for the nearest and next nearest neighbors, and treating the 
remainder of the array as point dipoles. The experimental H B shifts 
and the calculated interaction fields agree remarkably well (Fig. 4), 
confirming the validity of the mean-field description of the FORC- 
distribution and its quantitative predictive power, making the FORC 
technique a tool to extract numerical values of interaction fields. This 
is particularly important for disordered arrays where calculations of 
interactions are not easily achievable. 

Discussion 

In this work, systems of interacting nanomagnets were examined 
experimentally, numerically, and analytically, using the FORC 
technique. A mean-field analysis based on the concept of unmatched 
jumps accounted for all experimentally observed features of the 
FORC diagram, including its shifted ridge-and-edge structure and 
negative features. The tilting, shifting, and stretching of these struc- 
tures were identified as tools to extract quantitative information 
about the system, demonstrating the predictive power of the FORC 
technique. Construction of the FORC distribution through unmat- 
ched jumps, and recognizing the (de) magnetizing interactions as a 
particular case of (de) stabilizing interactions, presents an approach 
which can be used to evaluate any hysteretic system with the FORC 
technique. 

Methods 

Arrays of Co ellipses were fabricated by DC magnetron sputtering, in a vacuum 
chamber with a base pressure of 1 X 10" 8 Torr and Ar sputtering pressure of 2 X 10" 3 
Torr, on Si (100) substrates, in conjunction with electron beam lithography and lift- 
off techniques. Magnetic hysteresis loops were measured at room temperature using 
the magneto -optical Kerr effect (MOKE) magnetometer with a 632 nm HeNe laser 
having a 30 um spot-size, capturing the reversal behavior of —5,000 ellipses 32 . The 
magnetic field was applied parallel to the major axis of the ellipses. Each measurement 
was averaged over ~ 10 3 cycles at a rate of 1 1 Hz. The arrays were coated with a 60 nm 
ZnS layer to improve the signal-to-noise ratio 33 . FORC measurements were per- 
formed as follows 22 : from positive saturation the magnetic field is swept to a reversal 
field H R , where the magnetization M(H, H R ) is measured under increasing applied 
field H back to saturation, tracing out a FORC. The process is repeated for decreasing 
reversal field H R 30 . 

Ellipses were modeled as dipoles oriented parallel to their major axes. The inter- 
dipole spacing and magnetic moment per dipole in the 100 X 100 array were rep- 
resentative of the experimental system. Each dipole i was assigned an intrinsic 
coercivity H K * with a distribution experimentally determined from the sample having 
the weakest interactions, A3. The H in { dipolar interaction fields at dipole i were 
calculated on the mean-field level as aM(H), where a was calibrated such that aM s 
equals the analytically calculated H int at saturation. This mean-field formulation was 
extended by the first term of a cluster expansion, representing the nearest neighbor 
dipole interaction H nn ( explicitly: H in l = (xM(H)+H nn \ At each field step (AH = lOe) 
the total field H to l = H+H in l was compared to H K \ down-flipping occurred when 
H+H in t < —Hk and up-flipping occurred when H+H in t > H K \ until all dipoles 
became stable. 
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